


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1968 


Laminar incompressible flow in the entrance 
region of a square duct. 


Nielsen, David Joseph 


Monterey, California. Naval Postgraduate School 


http://ndl.handle.net/10945/12535 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
D U DLEY research materials and institutional publications created by the NPS community. 
get Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS'‘s first 
KNOX appointed — and published — scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 





NPS ARCHIVE 
1968 
NIELSEN, D. 


LAMINAR INCOMPRESSIBLE FLOW 
IN THE ENTRANCE REGION OF A SQUARE SUCT 


by 


David Joseph Nielsen 


TERRY 
Java POSTURSDUATE SAO 
waver, cus. 98040 











UNITED STATES 
NAVAL POSTGRADUATE SCHOOL 





THESIS 


LAMINAR INCOMPRESSIBLE FLOW 


IN THE ENTRANCE REGION OF A SQUARE DUCT 
by 


David Joseph Nielsen 


Marc h 196 3 








LAMINAR INCOMPRESSIBLE FLOW 


IN THE ENTRANCE REGION OF A SQUARE DUCT 


by 
David Joseph aba 
2 
Ensign, United States Navy 


B. S., United States Naval Academy 





Nid ee 


Submitted in partial fulfillment of the 
requirements for the degree of 
MASTER OF SCIENCE IN AERONAUTICAL ENGINEERING 
from the 
NAVAL POSTGRADUATE SCHOOL 


March 1968 





i Ge | 
ABSTRACT 
The development of the three dimensional, laminar velocity profile 

in the entrance length of a rectangular duct is investigated. The solution 
to this hydrodynamic problem is obtained from the full, incompressible 
Navier-Stokes equations and the continuity equation, in finite difference 
form, on the digital computer employing the computational method of 
Chorin (1). The solution yields the hydrodynamic velocities U, V, and 


W and the friction factor as a function of the distance from the entrance. 
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I, INTRODUCTION 

The design of compact heat exchangers for use in gas turbine 
regenerators, among other applications, entails the maximum uSe of 
the high heat transfer rates attainable in the entrance region. In this 
region the fluid is undergoing rapid transition from its uniform profile at 
the entrance to its fully established value farther downstream. A knowl- 
edge of the temperature profiles in the region is required to obtain the heat 
transfer rates. The temperature profile develops simultaneously with the 
velocity profile and is dependent upon it. The Graetz solution aSSumes a 
fully established velocity profile and thus eliminates most of the complexity 
of the hydrodynamic problem. This solution, however, fails to yield ac= 
curate results in problems involving large entrance region to overall length 
ratios. 

Various approximate methods have been devised for two dimensional 
channels, such as a circular pipe or parallel plate channels, to determine 
the development of velocity profiles. Until recently these solutions have 
invariably made boundary layer assumptions. Schlichting (5) obtained a 
solution for the flow between a pair of infinite parallel flat plates. He 
used two asymptotic series solutions, one based on Blasius’ solution of 
the boundary layer development expanded in the downstream direction, 
and the other based on the Hagen-Poiseuille solution of a parabolic 
velocity distribution, expanded in the upstream direction. He then 


joined these two solutions. 
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Recently Wang and Longwell (6) have employed the full momentum 
equations and the introduction of astream function to obtain a solution 
for the velocity development between parallel flat plates. The resulting 
equations were cast in finite difference form and the solution effected 
numerically. To date, no exact, three dimensional solutions for develop= 
ing velocity profiles have been obtained. 

The purpose of this analysis is to present an exact, finite- 
difference solution to the full incompressible Navier-Stokes equations 
for the hydrodynamic entrance region of a square duct. This solution is 
not nearly as simple as approximate methods. Its value, however, lies 
in its potential accuracy and the fact that it is an exact solution on 


which the credibility of various approximate solutions may be based. 
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II. ANALYSIS 


Be Governing Equations 


The governing equations are the momentum (Navier-Stokes) and the 
continuity equations. The complexity of the equations requires the fol- 
lowing assumptions concerning the flow: | 

1. steady and three dimensional 

2. laminar 

3. incompressible 
Through these assumptions, the energy and momentum equations become 
uncoupled, and may be expressed as follows: 

Case I: 


Momentum: 


du,us4 yy du 
Ye owe sm sii oP +7 [Ss + Zu. de) 











(1) 


M@iusv yy sv +wdy = 
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Be ud py WwW. = - 


| 
Tay ae eee (3) 


Continuity: 
QU ou + QW = 
3h " S94 s% 7 (4) 


1S 


Be Boundary Conditions 
The boundary conditions associated with the equations for a square 
duct with a uniform velocity at the entrance are: 
1, At the entrance (x = 0) 
v=w=0 
BS le 
2. At the walls (x=0, y= Youn 2a = Zea) 
u=v=w=0 
Cx Non-dimensionalization 


Equations (1), (2), (3) and (4) may be written in dimensionless form 


through the introduction of the following dimensionless variables: 





= Me! 
Use, Ve Wei Be oh, Ren SEE om 
yee Yoo Zot Tee 


The non-dimensional equations are then: 


Momentum: 


aX, IXY oe (6) 
V 
ay 7 3K> “SY soce (7) 


\Z SX* NY" Se (8) 


Continuity: 
gv , av, swe | 
WO - SY a iS =, (9) 
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D. Friction Factor 
To obtain an equation for the friction factor, the relation between 


the shear stress and the flow variables must first be introduced: 


a, Be a 
= <a Or = — 
Le 7" ayues 7 adit coaa (10) 
Since the duct is Symmetrical, either derivative, 2 or 9 , May be 


used in equation (10). By introducing a new variable ) , where Neyry , 


equation (10) may be written: 
Su 
= —— Lil 
T= p 3M \yeo (11) 
The friction factor, F, may be related to the Shear stress by the equation. 


f= ot 


oo (12) 


By substituting equation (11) into (12) the desired equation for the friction 


factor is obtained: 


(13) 





By defining a new dimensionless variable, = 4 , equation (13) may 


be written in dimensionless form 


> 


(14) 
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Ill. METHOD OF SOLUTION 


A. Technique 
The recently formulated technique of Chorin (1) was adopted for the 
solution of the governing equations. The principle of the method of solution 
suggested by Chorin lies in the introduction of an artificial compressibility, 
C , into the equations of motion in Such a way that the final results do not 
depend on ¢ . The solution is in essence a relaxation method. In employ- 


ing this method, it is necessary to introduce a modified set of governing 


equations: 
Momentum: 
2. you all ee a 
Rl By wal x x AY" OP (15) 


SY Oy Oven (16) 


IW Roly E 3W +4 a 264 Ts, Li . vW 


a) aX ay ee yee (17) 
Continuity: 

Ser er =o (18) 
State 

p= Os (19) 
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In this set of auxiliary equations © is an artificial density ,§ 
is an artificial compressibility, and P= C/s is an artificial — 
of state. 
Equations (15), (16), (17) and (18) may be put into finite difference 


form by introducing the following finite difference approximations: 


aU. UK KN, 2.7) UO AK, YET) Ub") = UK) 
3X AX 7 AY 


SU _ UCK+OX, Y, 2,7) HU 4K, Y,2,7) - U(KY, 2,7? AT) - 0 (KY,2, TAT) 
a 


_ vit) + ue) - UT )-ur) 
Ave 
The set of modified equations put into finite difference form and 


solved for the appropriate variables yields the following computing 


equations: 
US ee | RO 
) 13ST + 207 = Rey (X") =U ics) 


- RAT I yy)y(v) -ulY) viv] Bet Re AT] wea we) “Ue 


eee, u(x) +UK)- UT | +2 257) uty )+ UY) - UT | 


me 


+ BF [ueeue)-ury]-2E few)-Cor)]+ vr (20) 


NS, 


©) 


0)" aaa Sef )V(x*) = U(X vor) 
- Real | yi Y") Vw) Beal [i(e*)v(2) - We) WE’ ule") 


+a (x*) #90) - -v(r)}% aa F, | or ov (Y)-vit? 


+ 28h lve)evied-vr | sa ewv)-€(r)} er | (21) 


W(T)= 


ET RETPETS , net lu (x*)\W(x*) Ue WE)| 
~ Rell yy (y*)w(V*) -Wy rw(r)}- § Real |, 2") -W2 | 
+ 280 ate eae) -w0rg] #28 wren -WE 


+ ae we wie") +w(e) w(t v)/-35 ete)-ee)} wt 7 (22) 
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(23) 


B. Boundary Conditions 


In dimensionless form the boundary conditions for equations (15), 
(16), (17); and (18) are: 
1. At the entrance: (X = 0): 
V=W=0 
U=1 
2. Atthe walls Y=0,1 Z=0,1): 
W=V=U=0 
3. At X =X entrance length: 


W=V=0 


eS Veunly established 





Not 
_ OF Rees ae [1 SoSH) cos NZ 
io 2 eC) ” COsH(NT/2) (Nz) 


oes Stability . 


With the introduction of the modified equations (15), (16), (17), 


(18), and (19) an artificial Mach number is also introduced: 


Vo a ; Ya 
M = Re § (Seer an) 
Chorin (1) has stated that, in order to insure stability, the flow 


Mach number must be kept less than one. An additional requirement that 


must be met for stability of the set is that: 


Me cee CANN ov AY or AZ) - 


Zi 


D.. Calculation of F 
In order to calculate the friction factor, it is necessary to 


introduce a forward difference approximation as follows: 


i U(A+AR)-UMH) _ ulA)- UM) 








oan AN 
Using this approximation, equation (14) may be written: 
Spe URS 
e AY INFO 


Since equation (14) is applied at Damion U(Y\) is identically zero and 
is eliminated from the finite difference expression. 

Equation (24) may be used to obtain an average value of the friction ° 
factor. An average value of oF = may be determined at each X-wise 
location by summing U (H*) at each grid point on the boundary and 
dividing by the number of points used in the summation. The average 
value of the slope of the velocity profile may now be used to calculate 


an average friction factor as follows: 


F = Nam Re 2. Uwe (25) 
\) = number of points used 
oe = velocity at one grid point from the side wall boundary 


The grid system and summation points are shown in Figure |. 


Ee Numerical Procedure 
A finite difference solution was obtained, using the IBMsystem 360 


digital computer. Equations (20), (21), (22), and (23) involve approximating 


ZZ 


a solution for the entire flow field at time equal to T and T and then 

calculating a new approximation for time equal to jae The computed 

values of U, V, and W are tested for continuity in equation (23), and 

if eT) -€ al a does not equal zero, then a steady state solution 

has not yet been achieved. T 1s then set equal to T’ and T™ equal to T. 

Then the process is repeated until a steady solution is reached. That 

is until 7 = © - Lhis sequence serves to relax the set of equations. 
solution also requires that Mr be approximated at the grid points 

on the side wall boundary since equation (23) cannot be used at these 

points. To apply equation (23), one would have to impose the condition 

that the flow parameters are identical on either side of the channel walls. 

This assumption, however, is unrealistic. Instead, the following approxi- 


mation is introduced at the side wall boundary: 
CRY s — VVReeerr) 


Using the stability requirements, the solution to equations (20), (21), 

(22), and (23) was obtained using the following values of grid parameters: 
AX = .00164 AT =,.00000) AY= AZ=.1| Cama 

The value of AX was obtained using the entrance length of .0328, 

computed by McComas (4), for a hydraulic diameter of 1.0 anda 

Reynolds number of 1.0. The fully established velocity profile at this 

entrance length was taken from Knudsen (2). In this series solution for 

the fully established value of U, the pressure gradient Pry was replaced 

by e Or , where F is the fully established factor, also taken from 


Knudsen (2). 


The value of F computed by Knudsen (56.24) agrees well with the 
value suggested by Lundgren (3) (56.908) and McComas (4) (56.908). 
The computer program written in Fortran IV is included in the 


appendix. 
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WI. RESULTS 


Various combinations of AX,AY,A2,AT, and § were used in the 
solution, all of which satisfied the stability requirements set forth by 
Chorin (1). It was also found, however, that unless the ratio AT/, x? 
was of the same order of magnitude as the other ratios of this type, the 
solution would diverge. The values of AX,AY,Az, AT, and § given 
previously were found to be the most rapidly convergent values tried. 
Profiles were obtained using these values and Re equal to 1.0. They 
are plotted in Figure 2. The development of the centerline velocity is 
shown in Figure 3 and Table 1. Contrary to expectations, it was found 
that for any value of entrance length, the flow would not become fully 
established before it reached the x, imposed upon it. This was true 
even though the 2 was much larger than the value computed by McComas 
(.0328). In order to test the validity of the method, a solution for fully 
established flow in the duct was obtained using the fully established 
profile as boundary conditions at X = 0 and X = X entrance length. In 
this case, the profile at every X-wise location should be the same. 

Since the X-wise station X= Xe /2 is the last station to be effected 

by the relaxation technique, the profile at X= Xe/Q is plotted. For 200 
iterations the solution has not yet coriverged, and the probable error based 
on a calculation of U(Ke) /U(%$4) on the centerline is 15.4 per cent. After 
500 iterations, the probable error has decreased to 0.7 per cent, and after 
800 iterations. the probable error is reduced to 0.1 per cent. One can 


see that as the number of iterations is increagced, the accuracy of the 


NO 
1 


solution is increased until a very high degree of accuracy (99.93 per cent) 
is obtained at 1500 iterations. This is shown in Table II, Figure 4 and 
Figure 5. 

The average friction factor obtained as a function of X in the entrance 
length is shown in Figure 6 and Table III. The first four to six computed 
points were in error due to the finite value of AY and AZ _, which yield 
a finite value of shear stress at the entrance. In actuality the shear stress 
is infinite at the entrance, due to the discontinuity in the velocity at 
Vor Z=O. Inthe limit as AY approaches zero, equation (25) does, in 
fact, approach infinity. As the entrance length is approached, the value 
of F asymptotically approaches a value of 54. This value is in very good 
agreement with the fully established value that was introduced into the 
boundary conditions. It alSo is in good agreement with the value suggested 


by McComas (4) and Lundgren (3). 
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V. CONCLUSIONS 


As a result of this analysis, it may be concluded that: 

A solution using a finite difference approximation, and relaxation 
type of scheme suggested by Chorin (1), leads to a convergent and 
accurate solution albeit with substantially more computation than is 
inherent to an explicit scheme. The accuracy of the solution aera on 
the number of iterations made, the size of the finite-difference variables 
Ny Gan SY, AZ, and AT , and the value of the entrance length imposed 
upon the flow. This dependence of the solution upon the length of the 
entrance region is not what was expected and should be the subject of 
further investigation. The solution obtained in this analysis was based 
on a Reynolds number of 1.0. As one can see, the computing equations 
contain Reynolds number as a parameter and therefore depend upon it. 
Before a solution of this type may be used in the calculation of heat 
transfer in an entrance region, it would be necessary to investigate the 
behavior of this method at Reynolds numbers that are more typical of 
actual flow problems (R, = 10°) . A further limitation to the solution 
obtained in this analysis is that of computer time and storage space. 
This might be overcome by making use of the symmetry of the duct with 


a Suitable modification of the boundary conditions. 


if 
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APPENDIX A 


SAMPLE FORTRAN IV PROGRAM 
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TABLE II 
PROBABLE ERROR 
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€ AT X= 3 AT X= 
ITERATIONS Ding Pee bie X=X, ERROR (%) 
(COMPUTED) (ACTUAL) 
200 84575 1.0 15.4 
S ae 99288 1.0 ae 
800 99908 1.0 .09 
1! 00 99932 1.0 . 7 
(500 .99933 |. O pp 7 


39 


.00164 
.00328 
.00492 
-0065 6 
.00820 
-00984 
-OI1I148 
-O13S12 
.0 1640 
.0 1968 
.02296 
.0262 4 
.02788 
.02952 
.O3116 
.03280 


TABLE 111 
CALCULATED FRICTION FACTOR 


I RP QQ 


F 
8 0 
7T8.I| 
T 62% 
74.3 
72.6 
70.17 
69.085 
S (i> 
65.7 
64.16 
6i.1- 
586.7 
ro ro en 7 
5 4.82 
§ 4.64 
§ 4.4 
§ 4.2 


40 


INITIAL DISTRIBUTION LIST 


Defense Documentation Center 
Cameron Station 
Alexandria, Virginia 22314 


Library 
Naval Postgraduate School 
Monterey, California 93940 


Commander, Naval Air Systems Command 
Navy Department 
Washington, D. C. 20360 


Professor James A. Miller (Thesis Advisor) 
Department of Aeronautics 

Naval Postgraduate School 

Monterey, California 93940 


Chairman, Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 93940 


Professor A. E. Fuhs 
Department of Aeronautics 
Naval Postgraduate School 
Monterey, California 93940 


Dr. E. S. Lamar (Code 03C) 
Chief Scientist 

Naval Air Systems Command 
Navy Department 
Washington, D. C. 


Dr. kk. o. Bummgton 

Chief Mathematician 

Naval Air Systems Command 
Navy Department 
Washington, D. C. 


Commander 

Naval Ordnance Systems Command 
Navy Department 

Washington, D. C. 


41 


No. Copies 
20 


#0: 


Pe 


eZ 


ee 


14, 


Fo 


16. 


rae 


Office of Naval Research 1 
Navy Department 
Washington, D. C. 


Mr. G. L. Desmond l 
Aerodynamics and Structures Admin. (Code 320) 

Research and Technology 

Naval Air Systems Command 

Navy Department 

Washington, D. ©. 


Mr. I. Silver 

Propulsion Administrator (Code 330) 
Research and Technology 

Naval Air Systems Command 

Navy Department 

Washington, D. C. 


Office of Naval Research Ih 
(Attn: R. D. Cooper, Code 438) 

Mathematical Sciences Division 

Office of Naval Research 

Navy Department 

Washington, D. C. 


Office of Naval Research il 
Mathematical Sciences Division, Code 430 

Navy Department 

Washington, D. C. 


Office of Naval Research Il 
Air Programs Office 

Navy Department 

Washington, D. C. 


Ensign David J. Nielsen 3 
1510 N. Harmon 
Tacoma, Washington 98406 


Head, Engineering Department 1 


U. S. Naval Academy 
Annapolis, Maryland 21402 


4Z 


UNCLASSIFIED 


Security Classification 


DOCUMENT CONTROL DATA - R&D 


(Security claeelfication of title, body of ebetract and indexing annotation muet be entered when the overel! report ie claeeified) 


- ORIGINATING ACTIVITY (Corporate author) Za. REPORT SECURITY CLASSIFICATION 
Naval Postgraduate School 


Monterey, California 93940 
Laminar Incompressible Flow in the Entrance Region of a Square Duct 






















- REPORT TITLE 





- DESCRIPTIVE NOTES (Type of report end inclusive dates) 


. AUTHOR(S) (Last name, firet name, initial) 


NIELSEN, DAVID J., ENSIGN, USN 


6. REPORT DATE 7a. TOTAL NO. OF PAGES 


‘8a. CONTRACT OR GRANT NO. | 9a. ORIGINATOR'S REPORT NUMBER(S) 





BR PROJECT NO. 


9 b. QTH ER AER CRT NO(S) (Any other numbere thet may be aeeaigned 
le repo 


10. AVAILABILITY/LIMITATION NOTICES 


11. SUPPL EMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY 


19. ABSTRACT | 


The development of the three dimensional, laminar velocity profile in the 
entrance length of a rectangular duct is investigated. The solution to this 


hydrodynamic problem is obtained from the full, incompressible Navier-Stokes 
equations and the continuity equation, in finite difference form, on the digital 
computer employing the computational method of Chorin (1). The solution yields 
the hydrodynamic velocities U, V, and W and the friction factor as a function of 
the distance from the entrance. 





DD 52%. 1473 43 UNCLASSIFIED 
Security Classification 


UNCLASSIFIED 


Security Classification 


14 LINK A 


Laminar Flow 
Velocity Profile Development | 
Channel Flow 


Solution of Navier Stokes Equations 


‘J 
cd 1) » ‘ ® ? De 
SAY one * i. ss 1 fyeeee =O of o r " A on ee Oe 
% Codewtae vv — ccieiaie = - aided - os ere diate F] 
| # nolgpiete y ene tow ‘OR ; Sy MMT 
: ss Oa ages ee = et ie ME A - Ree —e Ee 
ak Ft 


DD now esd 473 (BACK) 44 UNCLASSIFIED 


Security Classification A-ipaes 























